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Abstract 

Statistical features of homogeneous, isotropic, two-dimensional turbulence is discussed on the 
basis of a set of direct numerical simulations up to the unprecedented resolution 32768^. By 
forcing the system at intermediate scales, narrow but clear inertial ranges develop both for the 
inverse and for direct cascades where the two Kolmogorov laws for structure functions are, for the 
first time, simultaneously observed. The inverse cascade spectrum is found to be consistent with 
Kolmogorov-Kraichnan prediction and is robust with respect the presence of an enstrophy flux. 
The direct cascade is found to be more sensible to finite size effects: the exponent of the spectrum 
has a correction with respect theoretical prediction which vanishes by increasing the resolution. 
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"Now, listen to me. You are living on a Plane. What you style Flatland is the 
vast level surface of what I may call a fluid on, or in, the top of which you and 
your countrymen move about, without rising above it or falling below it." 
Flatland by E.A. Abbott 

The existence of two quadratic inviscid invariants is the most distinguishirig feature of 
Navier Stokes equations in two dimensions. On this basis, R.H. Kraichnan 1[ predicted 
many years ago the double cascade scenario: when the turbulent flow is sustained by an 
external forcing acting on a typical scale if, an inverse cascade of kinetic energy E = l/2{v'^) 
to large scales {i ^ if) and a direct cascade of enstrophy Z = 1/2 (w^) to small scales 
{i ^ if) develop. In inverse and direct ranges of scales the theory predicts the kinetic 
energy spectrum E(k) ~ e'^^^k~^^^ and E(k) ~ rf/^k~^ with possible logarithmic corrections 
([2I). Here e and rj ~ /c^e are respectively the energy and the enstrophy injection rate. 

Navier-Stokes equations in two dimensions are now the prototypical model for turbu- 
lent systems displaying a double cascade scenario. From two-dimensional magneto-hydro- 
dynamics, to many geophysical model (such as Charney-Hasegawa-Mima), to wave tur- 
bulence models, the picture originally developed by Kraichnan has found many fruitful 
applications. 

Despite the expansion of the fields of applicability, it is remarkable that the verification of 
Kraichnan's theory, after more than 40 years from its formulation, is still partial. This is due 
to several reasons. First of all, the difficulties to generate a laboratory flow which is truly two 
dimensional on a large range of scales, limits the experimental approaches. From a numerical 
point of view, the situation in two dimensions is apparently very convenient with respect 
to three dimensions. A deeper analysis shows that this is not the the simultaneous 

simulation of two inertial ranges requires very large resolutions. Moreover, because time 
step is proportional to grid size, the computational effort for simulating two-dimensional 
turbulence can be even larger than in the three dimensional case. 

In the present paper we report numerical results on the statistics of the two cascades 
of two-dimensional turbulence on the basis of very high resolution (up to 32768^) direct 
numerical simulations. Together with previous results at lower resolutions (already reported 
on jsl) we obtain strong indications that the classical Kraichnan scenario is recovered in the 
limit of two inflnitely extended inertial ranges, although we are unable to address the issue 
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of possible logarithmic corrections in the direct cascade. 

The motion of an incompressible (V ■ u = 0) fluid in two dimensions is governed by the 
Navier-Stokes equations which are written for the scalar vorticity field a; = V x u as 

dtOJ + u ■ Voj = z/V^cu — au) + . (1) 

In ([1]) 1/ is the kinematic viscosity, is a forcing term and the friction term —auj removes 
energy at large scales in order to reach a stationary state. Alternatively, one can consider 
the quasi-stationary regime with a = in which the integral scale grows, according to 
Kolmogorov scaling, asL(t) = e^/^t^/^. In this case, Galilean invariant statistics (i.e. velocity 
structure functions or energy spectrum) is stationary at small scales I < L{t). We remark 
that the form of the friction term in ([1]) physically represents a crude approximation of the 



effects induced by bottom or air friction on a thin layer of fluid j4|. 

We numerically integrate ([I]) by means of a standard, fully dealiased, pseudo-spectral 
parallel code on a double periodic square domain of side = Ly = 2'k at spatial resolution 
up to = 32768. The forcing term in ([1]) is 5-correlated in time (in order to control 
energy and enstrophy input) and peaked on a characteristic forcing scale if. We use either 
a Gaussian forcing with correlation function (/^^(r, t)/tj(0, 0)) = F(5(t) exp(— (r/£/)^) or a 
forcing which has support on a narrow band of wavenumbers around kj = Tf/if in Fourier 
space. In both cases this ensures that energy and enstrophy input are localized in Fourier 
space and only a limited range of scales around the forcing scale is affected by the details of 



the forcing statistics. More complex 
direct effect on inertial range scales 



brcing, not localized in wavenumber space, can have a 



m 

. The forcing scale in all run is fixed at if = L^/ 100 
to allow the development of inertial ranges both at scales i > if (inverse cascade) and i < if 
(direct cascade). For the largest simulation run (A^ = 32768) we study the inverse cascade in 
the quasistationary regime with a = and we stop the integration of ([T]) when L{t) < Lx to 
avoid the pile- up of energy at the largest available scale. Table [T] reports the most important 
parameters for the simulations. 

The first information we get from the Table is related to the direction of the energy 
and enstrophy fluxes. According to the original idea of Kraichnan on the double cascade, 
in the ideal case of an infinite inertial range all the energy (enstrophy) injected should be 
transferred to large (small) scales. This can be thought as a limit case of a realistic situation 
in which the inertial range has a finite extension because of the presence of large and small 
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Label A B C D E 



N 2048 4096 8192 16384 32768 

V 2 X 10-5 5 X 10-6 2 X 10-6 1 X 10"^ 2.5 x 10"'^ 

a 0.015 0.024 0.025 0.03 0.0 

Ijjly 13 26 40 57 116 

Rx 7.9 15.4 21.5 26.0 36.0 

0.54 0.83 0.92 0.95 0.98 

r]ylr]i 0.96 0.92 0.90 0.88 0.98 

b 1.8 1.1 0.75 0.50 0.35 

TABLE I: Parameters of the simulations. N spatial resolution, v viscosity, a friction, Ij = n/kf 
forcing scale, Rx = Z'^^'^/rj^j Reynolds number for the direct cascade £i, = u^^'^/r]},^^ enstrophy 
dissipative scale, ej energy injection rate, friction energy dissipation rate (large-scale energy 
flux for run E), rfj enstrophy injection rate, iji, viscous enstrophy dissipation rate, 5 correction 
to the Kraichnan spectral exponent in the direct cascade. Viscosity is tuned to have for all runs 
kmax^u — 3. For the run E, for which a = 0, £a is kinetic energy growth rate. 

scale dissipation. The characteristic viscous scale and friction scale ia can be expressed 
in terms of the energy (enstrophy) viscous dissipation rate {t]^) and friction dissipation 
rate Ea (rja) by the relations il = e^/rh^a-i^d i"^ = Ea/rja- Energy and enstrophy balance 



equations in stationary conditions give 



BflB 



a 



rj^ ^ {Ui^f - 1 

r). 1 - {LItfY ^ ^ 

Therefore with an extended direct inertial range, iy ^ £/, one has Evjea 0, i.e. all the 
energy injected goes to large scales. Moreover, if ^ £j one obtains rja/rji, i.e. all the 
enstrophy goes to small scales to generate the direct cascade. Indeed, from Table [T] we see 
that increasing the resolution, i.e. Cf/Cy, the fraction of energy which flows to large scales 
increases. Because in our runs ^a/^f is constant with resolution and because £j, oc z/^/^ we 
expect that, according to (|2]), Sv/^a oc [it as indeed is shown in the inset of Fig. [T] 

Most of the enstrophy (around 90%) is dissipated by small scale viscosity. We observe a 
moderate increase of the large-scale contribution to enstrophy dissipation 7]^ by going from 



run A to D. This is a finite size effect because we have to increase the friction coefficient 
a with the resolution N in order to keep the friction scale ia constant when e„ 

grows. Indeed, for the run E without large-scale friction, the enstrophy flux to small scales 
almost balances the input. 










7^ 1 

' s 










-(c) >^ 




0.5 




A// 


to 

^> 
to 










10-2 












-1 


.(a) 






1 0'^ 1 Q-^ 

V 



1 10^ 10^ 10^ 10"^ 




1 10^ 10^ 10^ 10"^ 

k 



FIG. 1: (Color online), (a) Energy and (b) enstrophy fluxes in Fourier space for the runs of TableUl 
Fluxes for runs D and E are computed from a single snapshot. Inset (c): ratio of viscous over 
friction energy dissipation versus kinematic viscosity for the 5 runs, the line represents a linear fit. 

Figure [1] shows the energy and enstrophy fluxes in Fourier space defined as IlE{k) = 
—dt E{k')dk' and Iiz{k) = —dt k''^E{k')dk' (where E{k) is the energy spectrum and 
the time derivative keeps the contribution from nonlinear terms in ([T]) only Qj). We observe 
that, because resolution is changed by keeping ia ^ if constant, the only effect of increasing 
resolution on the inverse cascade is the growth of Ea/si (i.e. IlE{k)/ej) as discussed above, 
while the extension of the inertial range does not change. Despite the limited resolution of 
the inverse cascade inertial range {kf = 100), we observe an almost constant energy flux 
which develops independently on the presence of a direct cascade inertial range (run A). Of 
course, because of the presence of the two energy sinks (viscosity and friction) a plateau 
indicating a constant energy flux is clearly observable for the largest resolution simulation 
E only. On the contrary, the direct cascade does not develop for the small resolution runs 
as the dissipative scale is very close to the forcing scale (see Table [T]). A constant enstrophy 



flux Ilz{k) which extends over about one decade is on the other hand obtained for the most 
resolved run E. 

The behavior of the fluxes around k c:^ kf depends on the details of the injection: transi- 
tion from zero to negative (positive) energy (enstrophy) flux is sharp in the case of forcing 
on a narrow band of wavenumber (run D and E) while it is more smooth for the Gaussian 
forcing which is active on more scales. 
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FIG. 2: (Color online). Energy spectra for the simulation of Table |T] compensated with the inverse 
energy flux. Lines represent the two Kraichnan spectra Ck~^^^ (dashed) with C = 6 and 
(dotted). The inset shows the correction 5 to the Kraichnan exponent for the direct cascade 3 
obtained from the minimum of the local slope of the spectra in the range kj < k < as a function 
of the viscosity. Errorbars are obtained from the fluctuations of the local slope. The line has a 
slope 0.38 and is a guide for the eye. 



In Fig. |2] we plot the energy spectra of the different runs compensated with the energy 
flux. In the inverse range k < kj a, Kolmogorov spectrum E{k) = C6a^^k~^^^ is clearly 
observed for all simulations. The value of the Kolmogorov constant C ~ 6 is compatible 



with those obtained from more resolved inertial range [10| and it is found to be independent 
on the resolution. For what concerns the direct cascade, the spectrum is steeper than the 
Kraichnan prediction k~^. This effect is due to finite size effects, as it reduces by increasing 
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the resolution. In order to quantify the recovery of the Kraichnan exponent, we computed 
for all runs the local slope of the energy spectra in the range of wavenumber kj < k < ky. 
A plateau for the slope in this range of scales defines the scaling exponent —(3 + 5) of the 
energy spectrum in the direct cascade. In the inset of Fig. [2] we plot the measured value 
of the correction 5 as a function of the viscosity of the run. It is evident that, despite the 
fact the classical exponent —3 is not observed, the indication is that it should eventually be 
recovered in the infinite resolution limit z/ — t- 0. It is interesting to observe that for the most 
resolved run E, for which the enstrophy flux is almost constant over a decade of wavenumbers 
(see Fig. [T]), the exponent of the energy spectrum still has a significant correction 5 ~ 0.35. 
We remark that a clear observation of Kraichnan k~^ spectrum in simulations is obtained 



12|, while steeper spectra has also been 



using some kind of modified viscosity only 

m 

observed in simulations of ([1]) with a large scale forcing, i.e. resolving the direct cascade 
only [l3^. Therefore also for the direct cascade our simulations support the picture for which 
the statistics of one cascade is independent on the presence of the other cascade. 
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FIG. 3: (Color online). Third-order longitudinal velocity structure function S^{r) = ((5u||(r))^) 
for run E at final time. The two lines represent the Kolmogorov laws (H]) (continuous) and ([5]) 
(dotted). Inset: compensation of 5'3(r) with e^r (circles) and with Tj^r^ (triangles). Lines are the 
coefficient 3/2 (continuous) and 1/8 (dotted). 

We now consider small scale statistics in physical space, starting from velocity structure 
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function Sp{r) = {{Su\\{r)y) (with {Su\\{r) = (u(x + r) — u(x)) ■ r/r). The Kolmogorov 
relation, a consequence of constant energy (or enstrophy) flux in the inertial range, together 
with assumptions of homogeneity and isotropy, gives an exact prediction for the third-order 



longitudinal velocity structure function 6*3 (r) 14j-|l6|. For the inverse cascade it predicts 



3 

5*3 (r) = -Ear for r if (4) 

while for the direct cascade 

'^3('^) = \"nvr^ for r < (5) 
o 

The third-order velocity structure function for the simulation E is shown in Fig. |3l Both 
Kolmogorov laws are clearly visible with the predicted coefficients. We remark that this is 
the first time that the two fundamental laws (jl]) and (E]) are observed simultaneously. 

In Fig.HJwe plot velocity structure functions of different orders together with the compen- 
sation with Kolmogorov scaling 5'p(r) ^ (e^r)^/^. Although the range of scaling is very small, 
the presence of a plateau in the inverse cascade range of scales confirms that intermittency 
corrections are very small or absent in the inverse cascade range jioj ]. 

Velocity structure functions are trivially dominated by the IR contribution in the direct 
cascade range. Therefore to investigate higher order statistics of the direct cascade one 
has consider either increments of velocity derivatives (e.g. vorticity increments) or velocity 
second-differences (the latter having the advantage of being Reynolds-number-independent 
in the limit of zero viscosity). A Kraichnan energy spectrum would correspond dimen- 
.ionally to Eat vort.cty st™ctu.e function, and indeed .e.o =caUng exponents fo. p > 3 



[l7l | or logarithmic structure functions 18| are predicted in the limit of vanishing viscos- 
ity. Power-law intermittency corrections in the direct cascade of Navier-Stokes equations 
are excluded by theory (while logarithmic corrections are in principle possible), but it is 



known that the presence o 



generate intermittency 



19|-l21| 



inear friction terms in ([T]) can both steepen the spectrum and 



We have seen that in our simulations, even at highest resolution and without friction, we 
observe a correction to the spectral exponent, and therefore we cannot expect to observe 
theoretically predicted vorticity structure functions. Moreover, because of the limited res- 
olution of the direct cascade, no clear scaling in vorticity structure functions is observed. 
Nonetheless we can address the issue of intermittency by looking at the probability density 
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FIG. 4: (Color online). Longitudinal velocity structure function Sp{r) of order p = 2 (red squares), 
p = 4 (blue circles) and p = 5 (pink triangles) from run E and compensated with Kolmogorov- 
Kraiclman prediction (eQ.r)^/'^. 



functions of fluctuations of vorticity at different scales within the inertial range. The result, 
for run E is shown in Fig. |5] for both velocity and vorticity increments. For what concerns 
Su\\{r) we observe self-similar pdf in the inertial range of scales, in agreement with the nor- 
mal scaling of Fig. HI The shape of pdf is very close to Gaussian with a flatness which is 
around 3.3. On the contrary, vorticity increments 6uj{r) are definitely far from Gaussian 
with tails which are found to be very close to exponentials. Nonetheless, the shape of the 
pdf does not change substantially in the range of scale of the direct cascade, an indication 
of small intermittency also in this case. 

Velocity increments pdf in Fig. [5] cannot be exactly Gaussian as the energy flux, pro- 
portional to Ssi^r), requires a positive skewness. Energy and enstrophy fluxes are defined 



in physical space in terms of filtered fields, as described in |22|, |23| . We introduce a large 
scale vorticity field = * w and a large scale velocity field u,. = * u obtained from 
convolution with a Gaussian filter (^^.(x). From those fields, energy and enstrophy fluxes 
Ilr ' (x, t), representing the local transfer of energy/enstrophy from scales larger than r to 
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FIG. 5: (Color online). Probability density functions (pdf) of velocity longitudinal increments 
(5?i||(r) at scales r = 2.5if (red outer line) r = b.Oif (green middle line) and r = 10.0^/ (blue inner 
line) rescaled with rms values (left panel). Pdf of vorticity increments 6uj{r) at scales r = 0.2£f 
(red inner line) r = OAij (green middle line) and r = 0.8^/ (blue outer line) rescaled with rms 
values (right panel). Black curves are standard Gaussian. Data from run E. 

scales smaller to r, are defined as 

Ui''\^,t) = -{T^p)rVa{vp)r (6) 

ni^\^,t) = -{a^)rVaCOr (7) 
where iTaf3)r = (^^a^^/3)r " (^^a)r (^^/3)r and Mr = {VaUj)r " iVa)r(^r- 

Figure [6] shows the physical fluxes averaged over space at flnal time of simulation E. 
The two range of scales for the energy and enstrophy cascades are evident for r/£j > 1 and 
r/if < 1 respectively. The flnite mean values of fluxes are the results of strong cancellation: 
the ratio between the (absolute) mean value and the standard deviation at the scales ri ~ 
2.2if and r2 ~ 0.19£f corresponding to the peaks of the two fluxes are —0.19 and 0.16 
for energy and enstrophy respectively. The correlation among the two fluxes is small: the 
correlation coefficient between Ilrf ■* and Ilrf ■* is only C{ri, — —0.17 conflrming the picture 
of independence of the fluxes in physical spaces already observed at lower resolution ^]. 
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FIG. 6: (Color online). Average energy (lower line, red squares) and enstrophy (upper line, blue 
circles) fluxes in physical space for run E normalized with energy/enstrophy inputs. 

In conclusion, on the basis of very high resolution numerical simulations, we obtain strong 
evidence that the double cascade theory developed by Kraichnan more than 40 years ago 
is substantially correct. This result required massive resolution as two inertial ran ges have 



to be resolved simultaneously. It is worth remarking that, despite some effort 2^, |25|, the 
clear observation of the two cascade is still lacking in experiments. We hope that our results 
will stimulate further experimental investigations of the double cascade scenario. 

Numerical simulations has been performed within the DEISA Extreme Computing Ini- 
tiative program "Turbo2D". 
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